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Abstract 

We discuss several ways of illustrating fundamental concepts in statistical and thermal physics 
by considering various models and algorithms. We emphasize the importance of replacing students' 
incomplete mental images by models that are physically accurate. In some cases it is sufficient to 
discuss the results of an algorithm or the behavior of a model rather than having students write a 
program. 
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I. INTRODUCTION 



Mathematics is both the language of physics and a calculational tool. For example, the 
statements a = F/m and V • B = express the ideas that acceleration is the result of 
forces and magnetic field lines exist only as closed loops. The fundamental thermodynamic 
relation dS = {l/T)dU + {P/T)dV — {fi/T)dN implies that there exists an entropy function 
that depends on the internal energy U, the volume V, and the number of particles N. It 
also tells us that the temperature T, the pressure P, and the chemical potential /i determine 
how the entropy changes. 

Although textbooks and lectures describe the meaning of mathematical relations in 
physics, students frequently do not understand their meaning because there are few related 
activities that undergraduate students can do. Instead students frequently use mathemat- 
ICS clS db calculational tool for problems which lead to little physical understanding. To 
introductory students a = F/m is just one of many algebraic relations which needs to be 
manipulated. For more advanced students it is a differential equation which needs to be 
solved. 

The availability of symbolic manipulation software and calculators that can do much of 
the mathematical manipulations which students have been traditionally asked to do means 
that instructors need to think carefully about what are the appropriate activities for stu- 
dents. How many of the problems at the back of textbook chapters teach useful skills or 
help students learn physics? Which skills are important? How much, if any, understanding 
is lost if more of the calculations are done with the aid of a computer? Is most of the physics 
in the setup of the problem and in the analysis of the results? Should we spend more time 
teaching students to use software packages more effectively? 

In addition to challenging the way we teach the traditional physics curriculum, there is a 
growing interest in including more computation into the curriculum. Many physics teachers 
view computation as another tool analogous to mathematical tools such as those used to 
solve algebraic and differential equations. The theory is the same and computation is needed 
only when exact solutions are not available. Because the latter provide useful illustrations of 
the theory, computation need be only a minor part of the undergraduate (and even graduate) 
curriculum. This situation is the common practice at most institutions. 

A popular rationale for more computation in the curriculum is that it allows us to do 
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more realistic problems, is important in scientific research, and provides a useful tool for later 
employment. These reasons are all valid, but they might be less important than they seem. 
Although computation can allow us to consider more realistic problems, the consideration 
of such problems usually requires a much greater understanding of the system of interest 
than most undergraduates have the background and time to learn. Knowing a programming 
language is useful in employment, but might become less important as more and more work 
is done by higher level languages targeted toward specific applications. And even though 
computation is ubiquitous in research, more of it is being done using software packages. 
Just as few experimentalists need to know the details of how electronic instrumentation is 
constructed, few scientists need to know the details of how software packages are written. 
Moreover there usually is not enough time for students to write their own programs except 
in specialized courses in computational physics and simulation. 

In this paper we argue that computation should be incorporated into the curriculum 
because it can elucidate the physics. As for mathematics, computation is both a language 
and a tool. In analogy to models expressed in mathematical statements there are models 
expressed as algorithms. In many cases the algorithms are explicit implementations of the 
mathematics. For example, writing a differential equation is almost the same as writing a 
sequence of rate equations for the variables in a computer program. An advantage of the 
computational approach is that it is necessary to be explicit about which symbols represent 
variables and which represent initial conditions and parameters. 

In other cases the algorithm does not look like a traditional mathematical statement. For 
example, the Monte Carlo algorithms used in statistical physics do not look hke expressions 
for the partition function or the free energy. There also are models such as cellular automata 
that have no counterpart in traditional mathematics. 

In the following we will discuss some examples of how the consideration of algorithms 
and models can help students and instructors understand some fundamental concepts in 
physics. Because computation has had a great influence on statistical physics, we will focus 
our attention on this area. It is also the area in which we have the most expertise. 
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II. APPROACH TO EQUILIBRIUM 

A basic understanding of probability is necessary for understanding the statistical behav- 
ior underlying thermal systems. Consider the following question, which can be stated simply 
in terms of an algorithm.^ Imagine two bags and a large number of balls. All the balls are 
initially placed in bag one. Choose a ball at random and move it to the other bag. After 
many moves, how many balls are in each bag? Many students will say that the first bag will 
have more balls on the average.^ In this case it might be useful to show a simulation of the 
ball swapping process.^ It is then a good idea to ask students to sketch the number of balls 
in one bag as a function of time before showing the simulation. We can then discuss why the 
mean number of balls in each bag eventually becomes time independent. In particular, we 
can illustrate the presence of fluctuations during relaxation and in equilibrium and how the 
results depend on the number of particles. We can illustrate many of the basic features of 
thermal systems including the concepts of microstate and macrostate, the history indepen- 
dence of equilibrium (the equilibrium macrostate is independent of the initial conditions), 
ensemble and time averages, and the need for a statistical approach. We can also discuss 
why fluctuations in macroscopic systems will be negligible for most thermal quantities. 

The approach to equilibrium can be repeated with a molecular dynamics simulation (see 
Sec. Ill) in which particles are initially confined to one part of a container. After the release 
of an internal constraint the particles eventually are uniformly distributed on the average 
throughout the container. It is important in this example to show the students the basic 
algorithm and some simple computer code, so that they are convinced that there is no 
explicit "force" pushing the system toward equilibrium, but rather that equilibrium is a 
result of random events. The idea is to make explicit that the general behavior of systems 
with many particles is independent of the details of the inter-particle interactions. 

III. MOLECULAR DYNAMICS 

If students are asked what happens to the temperature when a gas is compressed, they will 
likely say that it increases. Their microscopic explanation will likely be that the molecules 
rub against each other and "give off heat."^'^ These students have a mental model which is 
similar to a system of marbles with inelastic coUisions. 
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This granular matter model of a gas and liquid is appropriate for the types of non- 
thermal macroscopic systems that students experience in everyday life such as cereals, rice, 
and sand. For molecular systems the coUisions are elastic, and a hard sphere model^ with 
elastic collisions is useful for understanding the properties of fluids and solids. Thus, student 
intuition has some value, but is limited in its ability to account for much of the phenomena 
of thermal systems. In particular, the only relevance of the temperature in hard sphere 
models is to set the the natural energy scale. 

A concrete and realistic model of thermal systems is provided by thinking about molecular 
dynamics with continuous inter-particle potentials for a system with a fixed number of 
particles and flxed volume.^ Each particle is subject to the force of every other particle in 
the system, with nearby particles usually providing most of the force. The dynamics of the 
ith particle of mass rrii is determined by Newton's second law, = Fnet,i/?7ii, which can be 
integrated numerically to obtain the position and velocity of each particle.^ Students should 
be asked to think about the appropriate inter-molecular force law and what would happen 
in the simulation and a real thermal system. For example, what is different about a gas 
and a liquid? This question will lead students to conclude that there must be an attractive 
contribution to the force law. Also systems do not completely collapse and thus there must 
be a repulsive part as well. Students can then be led to conclude that the force law must 
look something like the force law derivable from the Lennard-Jones potential. 

Next ask students about pressure, and you will likely not receive a coherent answer. Even 
if students know that pressure is force per unit area, it is likely that they will not be able 
to explain what that means in terms of a microscopic model of a gas. The simplest way 
to determine the pressure is to compute the momentum transfer across a surface. Without 
doing the simulation we can see that because momentum is proportional to velocity, faster 
particles will lead to greater momentum transfer per unit time. Then we can discuss what 
other macroscopic quantities increase with particle speed. Students will typically bring up 
temperature because they associate temperature with kinetic energy. From this discussion 
they can conclude that the pressure and temperature depend on the mean speed of the 
particles. Also, if the density increases at a given temperature, more particles will cross a 
given surface, and hence the pressure increases with the density. Students can then conclude 
that there are two independent quantities, the density and the temperature, which determine 
the pressure. 
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Because the total energy is conserved in an isolated container, students can understand 
that the total energy consists of both potential and kinetic energy. At this point some dis- 
cussion is needed to help students understand that there is potential energy that has nothing 
to do with gravity. Once that understanding is reached, it is easy to understand that the 
kinetic energy does not remain constant, and thus there will be small fluctuations in the 
temperature. Similar reasoning can lead to the idea of pressure fluctuations. Thus, by imag- 
ining the particles in a molecular dynamics simulation the concept of thermal fluctuations 
can be inferred. This result can be related to the fluctuations discussed in the ball swapping 
model in Sec. II. 

What is the role of temperature? We will discuss temperature again in Sec. IV, but here 
we consider how molecular dynamics can be used to think about its meaning. Imagine two 
solids such that initially the particles in one solid are moving much faster than the particles 
in the other solid. The two solids are placed in thermal contact so that the particles interact 
with each other across the boundary between the two solids. This model provides a concrete 
realization of thermal contact. Particles at the boundary will exchange momentum and 
energy with each other. The faster particles will give energy to the slower ones and energy 
will be transferred from the solid with the faster particles to the one with the slower particles. 
The net energy transfer will cease when all the particles have the same mean kinetic energy 
per particle, but not the same total energy per particle, potential energy per particle or mean 
speed. Thus by thinking about this system we can gain insight about the connection between 
kinetic energy and temperature. More importantly, we can emphasize that temperature is 
the quantity that becomes the same when two systems are placed in thermal contact. 

Our experience has been that a consideration of molecular dynamics in various contexts 
over several weeks leads to most students replacing their granular matter model by one in 
which energy is conserved. 

IV. AN IDEAL THERMOMETER 

In a previous paper^ we discussed the demon algorithm, which can be used as a model 
of an ideal thermometer and as a chemical potential meter. In a computer simulation the 
demon is an extra degree of freedom which can exchange energy or particles with a system. 
Energy is exchanged by making a random change in one or more degrees of freedom of the 
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system. If the change increases the energy, the change is accepted only if the demon has 
enough energy to give the system. If the trial move decreases the energy of the system, the 
extra energy is given to the demon. In a similar way the demon can exchange particles with 
the system. The energy distribution of the demon for a given number of particles in the 
system is the Boltzmann distribution from which the temperature can be extracted because 
the demon and system are in thermal equilibrium. If particle and energy exchanges are 
allowed, the distribution of energy and particles held by the demon is a Gibbs distribution 
from which the chemical potential can be extracted. 

A discussion of the demon is a good way of introducing the concepts of system, heat bath, 
and ensembles. The demon algorithm simulates the microcanonical ensemble because the 
total energy of the system plus the demon is fixed, and the demon is a negligible perturbation 
of the system. An alternative interpretation is that the demon is the system of interest and 
the remaining particles play the role of the heat bath as in a canonical ensemble. The demon 
is unusual because its energy is both the energy of the system and the energy of a microstate. 
Similarly, the state of the demon is both a microstate and a macrostate. In realistic thermal 
systems there are many microstates corresponding to a given macrostate, and thus there 
are distinctions between a macrostate, a microstate, and a single particle state. The demon 
algorithm provides a concrete example for discussing such concepts. 

V. MARKOV CHAINS AND THE METROPOLIS ALGORITHM 

One of the most common algorithms used to simulate thermal systems is the Metropolis 
algorithm. -"^^ In this section we discuss the underlying theory behind this algorithm, which 
will allow us to introduce concepts such as probability distributions, the Boltzmann distri- 
bution, Markov chains, the partition function, sampling, and detailed balance. In Sections 
VI- VIII we will extend these ideas to newer algorithms that are currently of much interest 
in research. 

Most students are familiar with probabilities in the context of dice and other simple games 
of chance, but they have difficulty understanding what probabilities mean in the context of 

thermal systems. 

Students are familiar with systems that evolve in time according to deterministic equa- 
tions such as Newton's second law. We can also consider stochastic processes such as those 
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used in Monte Carlo simulations (see Sec. II). In this case we consider an ensemble of 
identical copies of a system. These copies have the same macroscopic parameters, interac- 
tions, and constraints, but the microscopic states (such as the positions and momenta of 
the particles) are different in general. We imagine a probabilistic process which changes the 
microscopic states of the members of the ensemble. An example is a Markov process for 
which the probability distribution of the states of the ensemble at time t depends only on 
the probability distribution of the ensemble at the previous time; that is, the system has no 
long term memory. A Markov process can be represented by the equation: 

P{t + At)^7rP{t), (1) 

where P{t) is a column matrix of n entries, each one representing the probability of a 
microstate; tt is a matrix of transition probabilities. The ith entry of P{t) gives the fraction 
of the members of the ensemble that are in the ith state at time t. By repeatedly acting on 
P with the transition matrix we arrive at the distribution at any later time t. If a system is 
in a stationary state, then 

P{t + At)=P{t). (2) 

How can we be sure that the Markov process we use in our simulations will lead to the 
theoretically correct equilibrium distribution? We begin with the detailed balance condition 

where TCi^j is the probability of the system going from state i to state j, and Pi is the 
probability that the system is in state i. Detailed balance puts a constraint on the transition 
probabilities. As we will show, if the stochastic process satisfies detailed balance, then a 
system which is in a stationary state will remain in that state. 

We now show that this property is satisfied for an equilibrium system at temperature 
T = l/k/3 for which the probability of a microstate is given by the Boltzmann distribution. 
If we substitute the Boltzmann distribution Pi = e~^^^ /Z into Eq. (3), we obtain 

7r,_,-e-^^^=7r,-_,e-^^^ (4) 

where Ei is the energy of state i and k is Boltzmann's constant. Note that the partition 
function Z — Ylii^'^^^ does not appear in Eq. (4). This absence is important because Z is 
usually not known. 



We use Eq. (4) and the Boltzmann distribution for the desired stationary state to calculate 
the right-hand side of each row of Eq. (1): 

P,(t + Ai) = 5;]7r,_.,P,(i). (5) 

j 

We next use the detailed balance condition in Eq. (4) to obtain 

m + At) = ^ [7r,_,e-^(^-^^)] P,{t). (6) 
j 

We want to see if an ensemble that is described by the Boltzmann distribution will remain 
in this distribution. We replace Pj in Eq. (6) by the Boltzmann distribution e"^^^ /Z so that 

P.{t + At)^Yl [^^-^Je-^^^^-^^^ (7) 
j 

We can simphfy Eq. (7) by writing 

P,(i + Ai) = ^ Jj7r,_.,-, (8) 

j 

where we have taken the factor that does not depend on j out of the sum. Because the 
system has probability unity of going from the state i to all possible states j, the sum over 
j must equal unity, and thus Eq. (8) becomes 

P,(t + At)^^. (9) 

Thus, we have shown that if we start with an ensemble distributed according to the Boltz- 
mann distribution and use a transition probability that satisfies detailed balance, then the 
resulting ensemble remains in the Boltzmann distribution. Note that detailed balance is suf- 
ficient, but we have not shown that detailed balance is necessary for an ensemble to remain 
in a stationary state, and it turns out that it is not. 

Although detailed balance specifies the ratio ni^j/nj^i, it does not specify vTj^j itself. 
There is much freedom in choosing iTi^j, and the optimum choice depends on the nature 
of the system being simulated. The earhest and still a popular choice is tti^j — Ai^jWij, 
where Wij is the probability of making an arbitrary trial move and Ai^j is the Metropolis 
acceptance probability given by 

Ai^j = min |l, e-^(-^^-^*)}. (10) 



Equation (10) can be shown to satisfy detailed balance if = Wji (see Problem lb). For 
example, for an Ising model with N spins choosing a spin at random implies Wij — Wji — 
1/N. 

Our discussion of probability has been a mixture of some abstract ideas from probability 
theory and the example of the Metropolis algorithm. The discussion can be made more 
concrete by considering a specific model, such as the Ising model of magnetism or a Lennard- 
Jones system of particles. 



VI. DIRECT ESTIMATION OF THE DENSITY OF STATES 

The density of states g{E) is defined^^ so that g{E)AE is the number of microstates 
with energy between E and E + AE. For most students the density of states is an abstract 
quantity and many confuse the density of states of a many body system with the single 
particle density of states. The following discussion provides an algorithm which makes the 
meaning of g{E) more concrete. 

If the density of states g{E) is known, we can calculate the mean energy and other 
thermodynamic quantities at any temperature from the relation 

<e> = J2ep{e), (11) 

E 

where the probability P{E) that a system in equilibrium with a heat bath at temperature 
T has energy E is given by 

^ (12) 

and Z — '^e9(^) ^"^^ ■ Hence, the density of states is of much interest. 

Extracting the density of states from a molecular dynamics simulation or from the 
Metropolis and similar algorithms is very difficult. For example, suppose that we try to 
determine g{E) by doing a random walk in energy space by flipping the spins in the Ising 
model at random and accepting all microstates that are obtained. The histogram of the 
energy, H{E), the number of visits to each possible energy E of the system, would converge 
to g{E) if the walk visited all possible configurations. In practice, it would be impossible to 
realize such a long random walk given the extremely large number of possible configurations 
in even a small system. 
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Recently a Monte Carlo algorithm due to Wang and Landau^^ has been developed for de- 
termining g{E) directly. The idea is to simulate a system by making changes at random, but 
to sample energies that are more probable less often so that H{E) becomes approximately 
independent of E. The acceptance criteria is given by 

where g{E) is the running estimate of g{E). This acceptance criteria favors energies which 
have a low density of states, and because these states are visited more often, the result is that 
a histogram of the energy of the visited configurations will be approximately independent of 
E or "fiat." It is possible to show (see Problem Ic) that the probability of the ith microstate 
is proportional to l/g{Ei) in the stationary state. 

How do we implement Eq. (13) when we don't know g{E), which is the goal of the 
simulation? The second part of the algorithm is to make an initial guess for g{E) and 
then improve the estimate of g{E) as the simulation proceeds. The simplest guess is to let 
g{E) — 1 for all E and then to update g{E) and H{E) after each trial move: 

~gt+,{E) = rgt{E), (14a) 

or 

\ngt+i{E)=\ngt{E)+\nf, (14b) 

with 

Ht+,{E)^HtiE) + l, (15) 

and f > 1; t represents the number of updates. (The value of E in Eqs. (14) and (15) is 
unchanged if the trial move is not accepted.) Because g{E) increases rapidly with E, we 
need to use In ^(£') to implement this algorithm on a computer. 

The combination of Eqs. (13) and (14) forces the system to visit less explored energy 
regions due to the bias in the acceptance probability in Eq. (13). For example, if the 
current estimate g{E) is too low, moves to states with a lower value of g{E) have a greater 
probability. In this way the values oi g{E) will gradually approach the "true" g{E). 

Initially we choose / to be large (typically / = e ~ 2.7) so that g{E) will change rapidly. 
After many moves the histogram H{E) will become approximately flat, and we then decrease 
/ so that g{E) doesn't change as much. The usual procedure is to let / — > -^/J and continue 
the updates of g{E) and the changes in / until / — 1 ~ O(10~*). 
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As can be seen in Eq. (11) the only quantity that differentiates the mean energy and the 
specific heat of one system from another is the density of states. However, the density of 
states of many systems is quahtatively similar. This similarity is illustrated in Fig. 1, which 
shows the density of states for 64 Ising spins in one, two, and three dimensions. More insight 
can be gained from the plot in Fig. 2 of the probability P{E) superimposed on a plot of 
In g{E). The plot shows the range of values of E that are important for a given temperature. 
The competition between the increase of g{E) with E and the decrease of the Boltzmann 
factor leads to the peak in P{E). 

One of the interesting features of the Wang-Landau algorithm is that it samples states 
that are of little interest for thermal systems, unlike the Metropolis algorithm which rarely 
samples such states. For example, an application of the Wang-Landau algorithm to the Ising 
model would lead to a parabolic- like curve for \a.g{E) with a maximum at £^ = as shown 
in Fig. 1. Positive energy states, with most of the nearest neighbor spins of opposite sign, 
are present. Because large positive energy states would not be observed in a Metropolis 
simulation of the Ising model, it is clear that we are measuring a temperature independent 
property that depends only on the nature of the system. 

VII. DIRECT ESTIMATION OF T(E) IN THE MICROCANONICAL ENSEMBLE 

Another way to determine g{E) is to exploit the correspondence between the density 
of states and the thermodynamic temperature T{E) and to update the latter rather than 
g{E)}'^ The approach is based on the relation between the microcanonical entropy, S{E) — 
\ng{E) (with Boltzmann's constant k = 1), and the inverse temperature, p{E) = l/T{E) = 
dS/dE. The main virtue of this method of determining g{E) in the present context is that 
it is an application of the relation between the entropy and the temperature. 

For simplicity, we discuss the algorithm for the Ising model for which the values of the 
energy are discrete. We label the possible energies of the system by the index j with j — 1 
corresponding to the ground state and write the temperature estimate T{Ej) as Tj. After a 
trial move to state j we update the entropy Sj as for the Wang-Landau algorithm and also 
update the temperature estimate at j ± 1. From Eq. (14a) or 

Sj^t+i = Sj,t + ^f, (16) 
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we can write the central difference approximation for the inverse temperature (3 as 

= ^ [Sj+^, - Sj.^,]/{2AE), (17) 



dS 
dE 



E=Ei 



where AE is the energy difference between state j and states J ± 1. (For the Ising model 
on a square lattice, l\E — 4J. In general, we would choose a larger bin size so that several 
states would correspond to one bin.) The energy Ej in Eq. (17) is the value of the energy 
of the system after a trial move. On a visit to Ej we use the updated value of Sj and the 
unchanged values of Sj±2 to determine new estimates for the temperature. Prom Eq. (16) we 
hme j3j+i^t+i = [Sj+2,t-Sj^t+i]/i'2AE) =j3j^^-5f\ and^-.i^t+i = [Sj^t+i-Sj-2,t]/{'^AE) = 
(3j-i,t + 5f with 5f = ln//2A£^. Hence, we obtain 

^j±i,t+i = V/^iii = (^j±i,tTj±i^f, (18) 

where oij±i^t = T 5/^j±i,t]- Note that Tj_i^t+i is decreased and Tj^i^t-i is increased 
while Tj is unchanged. In this way T{E) will converge to a monotonically increasing function 
of the energy. It is best to restrict the updates to a finite range of temperatures between 

Tmin and Tmax- 

The acceptance probability of the trial moves is given by Eq. (13) so we also need to 
update the entropy as is done for the Wang- Landau algorithm. The values of / are changed 
as for the Wang-Landau algorithm. An example of the converged T{E) for the Ising model 
is given in Fig. 3. 

Once T{E) converges, the entropy estimate is found by integrating T{E) with respect to 
E. In the simplest interpolation method the entropy is given by 

S{E) = PjAE. (19) 

3 Jmin 

These values of S{E) are then used to obtain g{E), and the thermodynamic properties of the 
system are determined from Eq. (11). Because the continuum entropy S{E) can be obtained 
by integrating the interpolated Tj, the updating of T{E) as in Eq. (18) is especially useful for 
systems where the energy changes continuously. Even more interesting is the use of Eq. (18) 
with molecular dynamics simulations at various temperatures. ^'''^^ 
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VIII. DIRECT MEASUREMENT OF THE PARTITION FUNCTION 



Another Monte Carlo method combines the Metropohs algorithm with the Wang-Landau 
method to directly compute the partition function Z at all temperatures of interest. The 
method uses two types of trial moves: a standard Metropohs Monte Carlo move such as 
a flip of a spin which changes the energy at fixed temperature, and a move to change the 
temperature at fixed energy. The acceptance rule is given by 

where Ei is the energy of the current configuration, /9j is the current inverse temperature 
and is the current estimate of the partition function for inverse temperature For 

an energy move, /3j = and for a temperature move = Ej. At each trial move a decision 
to choose a temperature instead of an energy move is made with a fixed probability such 
as 0.1%. To update the partition function we use the same procedure as for the density of 
states in the Wang-Landau method: 

lnZ(/3) ^lnZ(/3)+ln/. (21) 

This procedure will lead to a stationary state with probability 

Pi = e^'^VZ(A). (22) 

A plot of In [Z{T)/Z{T = 0.1]/A^ for the two-dimensional Ising model with = 32 x 32 
spins is shown in Fig. 4. We rarely show plots of partition functions in a thermal physics 
course. The plot shows that In Z does not change very much at low temperatures, increases 
rapidly near the phase transition, and then increases slowly for higher temperatures. 

IX. SUMMARY 

Our focus has been on understanding important concepts in statistical mechanics by 
considering simulations of concrete models of thermal and statistical systems. Instructors 
can choose how to use simulations in their courses. For example, students can be asked 
to write programs with the use of templates. Another strategy is to have students run 
existing programs and modify some of the parameters and explain the results.^ 
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Computers are omnipresent in physics research. Much of this use is for data analysis, sym- 
bohc algebra (for example, to calculate Feynman diagrams), and the control of experiments. 
Statistical physics is an area where the development of new algorithms and simulations has 

qualitatively changed the kinds of systems and problems that can be considered. These 
developments make it even more important to think about the ways that computers should 
change the way we teach thermal and statistical physics. 

X. SUGGESTIONS FOR FURTHER STUDY 

1. Detailed balance. 

(a) We showed that the detailed balance condition in Eq. (4) ensures that the Boltz- 
mann probability distribution is stationary. Show that the detailed balance con- 
dition, Eq. (3), ensures that the distribution Pi is stationary. 

(b) Show that the Metropolis algorithm satisfies detailed balance if Wij — Wji. Show 
that the symmetric acceptance probability, 

also satisfies detailed balance. 

(c) Show that the stationary probability distribution for the Wang-Landau algorithm, 
Eq. (13), is Pi — c/g{Ei), where c is a constant independent of E. 

2. Approach to equilibrium. 

(a) Write a program that implements the balls and bags example discussed in Sec. II 
or use the application/applet available from (stp. clarku.edu). ^ 

(b) Plot the number of particles on the left side of the container as a function of time. 
How long does it take for a system of N — 64 balls to come to equilibrium? Then 
estimate this time for N — 128, 256, and 512 balls. 

(c) Once the system reaches equihbrium there are still fluctuations. Find a relation 
for the magnitude of the fluctuations as a function of N. 
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(d) Does there exist an initial configuration of the balls that comes to a different 
equilibrium state than the ones you have simulated so far? Why is the answer to 
this question important in statistical mechanics? 

3. Wang- Landau algorithm 

(a) Calculate by hand the density of states ior N — 6 Ising spins in one dimension 
with toroidal boundary conditions. The minimum energy is —6 and the maximum 
energy is +6. There are 2^ — 64 microstates. 

(b) Write a program that uses the Wang-Landau algorithm to determine the density 
of states for the one-dimensional Ising chain and compare your results with your 
hand calculation. 

4. Additional problems related to the topics in this article as well as many other topics 
in statistical physics can be found at (stp.clarku.edu) or EPAPS.^ 
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FIG. 1: Semi- log plot of the exact densities of states for AT = 64 spins for one, two, and three 
dimensions. The results for three dimensions were generated using the method discussed in Ref. 18. 
The results of the Wang-Landau algorithm are indistinguishable from the exact results for these 
small systems. Note the tiny deviations from a smooth curve at \E\ 120 for d = 2 and \E\ ^ 170 
for d = 3. These are signatures of a phase transition in the thermodynamic limit. For d = 1 the 
Ising model does not have a phase transition. 
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FIG. 2: The probability P{E) of the energy E for the Ising model on a 32 x 32 square lattice for 
various temperatures. Superimposed on the same plot is the density of states. 
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FIG. 3: The estimated energy-dependence of the temperature for the Ising model on a 32 x 32 
square lattice as determined by the method of Ref. 14. The simulation was done with Tmax = 4.5 
and Tmin = 1-3 using AE = 8 and converged to / = 1 -|- 10~^^ after 4 x 10^ Monte Carlo iterations 
per spin with the initial modification factor /o = 1.00005. The data was generated by Jaegil Kim. 
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FIG. 4: Simulation results for [\nZ{T)/Z{Q.l)]/N for a 32 x 32 Ising lattice. T = 0.1 is the lowest 
temperature simulated. The algorithm used and the temperatures simulated are the same as in 
Ref. 16. The final value of In / 3.16 x 10~^. The number of MC steps at each value of / was 
lOOAT/ In /. 
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